Predict with sims from density and diet models to propagate uncertainty
# Predict and then for loop to summarise and avoid having too large pred grids in the memorynsim <-500sim_dens <-predict(mcod, newdata = pred_grid, nsim = nsim) |>as.data.frame()sim_her <-predict(mher, newdata = pred_grid, nsim = nsim) |>as.data.frame() |>bind_cols(pred_grid |> dplyr::select(X, Y, year))sim_spr <-predict(mspr, newdata = pred_grid, nsim = nsim) |>as.data.frame() |>bind_cols(pred_grid |> dplyr::select(X, Y, year))sim_sad <-predict(msad, newdata = pred_grid, nsim = nsim) |>as.data.frame() |>bind_cols(pred_grid |> dplyr::select(X, Y, year))
The density prediction is used later for spatial overlap + it gets left_joined to the other sims, so it gets a special treatment
pred_dens_space <- pred_grid_density |>filter(year %in%c(1994, 2022)) |>summarise(density_mean =mean(sim_density),.by =c(year, X, Y) )annotations <-data.frame(xpos =c(-Inf, -Inf),ypos =c(-Inf, Inf),year =c(1994, 2022),hjust =c(-0.3, -0.3),vjust =c(-29.2, 1.5))# Max density by species (because I trim the plot)pred_dens_space |>summarise(max =max(density_mean))
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
ggsave(paste0(home, "/figures/spatial_cod.pdf"), width =17, height =9, units ="cm")
Plot annual predation with uncertainty
# Summarise predator biomass by year and sim and join that to calculate per capita predationpred_dens <- pred_grid_density |>summarise(cod_biomass =sum(sim_density *9), # area is 9 km2.by =c(year, sim) )# Get the predation by year in a loop, else it becomes too long with all the sims, running out of memory...herlist <-list()sadlist <-list()sprlist <-list()tic()for (i inunique(pred_grid$year)) { sim_her_sub <- sim_her |>filter(year == i) sim_sad_sub <- sim_sad |>filter(year == i) sim_spr_sub <- sim_spr |>filter(year == i) pred_grid_density_sub <- pred_grid_density |>filter(year == i) herlist[[i]] <-get_annual_predation(sim_her_sub, prey_name ="Herring") sadlist[[i]] <-get_annual_predation(sim_sad_sub, prey_name ="Saduria") sprlist[[i]] <-get_annual_predation(sim_spr_sub, prey_name ="Sprat")gc()}tot_pred <-bind_rows(bind_rows(herlist),bind_rows(sadlist),bind_rows(sprlist))toc()
269.777 sec elapsed
# Apply functions and combine# This is the old way when we did it with all years together... then we're limited by how long the df can be# tot_pred <- bind_rows(# get_annual_predation(sim_her, prey_name = "Herring"),# get_annual_predation(sim_sad, prey_name = "Saduria"),# get_annual_predation(sim_spr, prey_name = "Sprat")# )# Prepare data for fitting gams and for plottingclean_pred <- tot_pred |>mutate(species =as.factor(species))# Fit gamma-GAMs to annual predation metricsgamma_gam_pop <-sdmTMB(pred_median ~ species +s(year, by = species),spatial ="off",family =Gamma(link ="log"),data = clean_pred |>mutate(pred_median = pred_median /1000))gamma_gam_cap <-sdmTMB(cap_median ~ species +s(year, by = species),spatial ="off",family =Gamma(link ="log"),data = clean_pred)nd <-expand_grid(species =as.factor(c("Sprat", "Herring", "Saduria")),year =unique(clean_pred$year))p_pop <-predict(gamma_gam_pop, newdata = nd, se_fit =TRUE)p_cap <-predict(gamma_gam_cap, newdata = nd, se_fit =TRUE)p <-bind_rows( p_pop |>mutate(metric ="Population-level (tonnes)"), p_cap |>mutate(metric ="Per capita (kg/kg)"))pop <- clean_pred |> dplyr::select(year, species, pred_median, pred_lwr, pred_upr) |>rename(median = pred_median,lwr = pred_lwr,upr = pred_upr ) |>mutate(median = median /1000,lwr = lwr /1000,upr = upr /1000 ) |>mutate(metric ="Population-level (tonnes)")cap <- clean_pred |> dplyr::select(year, species, cap_median, cap_lwr, cap_upr) |>rename(median = cap_median,lwr = cap_lwr,upr = cap_upr ) |>mutate(metric ="Per capita (kg/kg)")clean_plot <-bind_rows(pop, cap)p0 <-ggplot(clean_plot, aes(year, median)) +facet_grid2(metric ~ species,switch ="y",scales ="free_y", independent ="y" ) +geom_errorbar(aes(ymin = lwr, ymax = upr),width =0, color ="grey40", alpha =0.7,linewidth =0.4 ) +geom_point(color ="grey30", alpha =0.8) +geom_line(data = p, aes(year, exp(est)), color = pal[2], linewidth =0.8) +geom_ribbon(data = p, aes(x = year,ymin =exp(est -1.96* est_se),ymax =exp(est +1.96* est_se) ),fill = pal[1], alpha =0.3, inherit.aes =FALSE ) +theme(aspect.ratio =6/7,strip.text.x.top =element_text(size =8.2, margin =unit(rep(0.06, 4), "cm")),strip.text.y.left =element_text(size =8.2, margin =unit(rep(0.06, 4), "cm")) ) +labs(x ="Year", y ="Relative predation")tag_facet(p0, fontface =1, col ="gray30", size =3)
ggsave(paste0(home, "/figures/pred_yearly.pdf"), width =17, height =9, units ="cm")
Plot weighted predation in space
pred_grid_density_sub <- pred_grid_density |>filter(year %in%c(1994, 2022))tic()spatial_pred <-bind_rows(get_spatial_predation(sim_her, years =c(1994, 2022), prey_name ="Herring"),get_spatial_predation(sim_sad, years =c(1994, 2022), prey_name ="Saduria"),get_spatial_predation(sim_spr, years =c(1994, 2022), prey_name ="Sprat"))toc()
17.11 sec elapsed
plot_map2 <- plot_map +guides(fill =guide_colorbar(position ="inside")) +theme(legend.position.inside =c(0.14, 0.84),legend.key.width =unit(0.2, "cm"),legend.key.height =unit(0.25, "cm"),legend.text =element_text(size =7),legend.direction ="vertical",legend.title =element_text(size =9),plot.title =element_text(margin =margin(0, 0, -10, 0)) ) +facet_wrap(~year, ncol =1)annotations <-data.frame(xpos =c(-Inf, -Inf),ypos =c(-Inf, Inf),year =c(1994, 2022),hjust =c(-0.3, -0.3),vjust =c(-20, 1.5))# Max predation by species (in case trim the plot)spatial_pred |>summarise(max =max(pred_mean), .by = species)
# A tibble: 3 × 2
species max
<chr> <dbl>
1 Herring 4.41
2 Saduria 84.7
3 Sprat 63.4
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
tag_facet(p1, fontface =1, col ="gray30")
ggsave(paste0(home, "/figures/supp/cv_density_relative_prey.pdf"), width =19, height =10, units ="cm")
Plot relative prey weight predictions in space for two years
# Here we can use the get_cv_predation and just use the median# Max predation by species (if I trim the plot)cv_pred |>summarise(max =max(mean), .by = species)
# A tibble: 3 × 2
species max
<chr> <dbl>
1 Herring 0.00649
2 Saduria 0.0325
3 Sprat 0.0324
# Plot!viridis(n =3, option ="magma") # high color